
#include <cstdio>
#include <cstring>
#include <cmath>
#include <cstdlib>
#include <iostream>
#include "lbfgsb.h"

using namespace std;

//const int SIXTY = 60;
//const int nmax = 40480, mmax = 20;

int main() {

    //c     This simple driver demonstrates how to call the L-BFGS-B code to
    //c       solve a sample problem (the extended Rosenbrock function 
    //c       subject to bounds on the variables). The dimension n of this
    //c       problem is variable.
	int nmax = 40480, mmax = 20;
    int SIXTY = 60;
    //int nmax = 1024, mmax = 17;
    //c        nmax is the dimension of the largest problem to be solved.
    //c        mmax is the maximum number of limited memory corrections.

    //c     Declare the variables needed by the code.
    //c       A description of all these variables is given at the end of 
    //c       the driver.
    //char task[SIXTY], csave[SIXTY];
	char* task = new char[SIXTY];
	char* csave = new char[SIXTY];

    long int lsave[4];
	long int n, m, iprint, isave[44];
	double f, factr, pgtol, dsave[29];
    //long int n, m, iprint, nbd[nmax], iwa[3 * nmax], isave[44];
    //double f, factr, pgtol, x[nmax], l[nmax], u[nmax], g[nmax], dsave[29];
            //wa[2 * mmax * nmax + 4 * nmax + 12 * mmax * mmax + 12 * mmax];
	long int* nbd = new long int[nmax];
	long int* iwa = new long int[3 * nmax];
	double* x = new double[nmax];
	double* l = new double[nmax];
	double* u = new double[nmax];
	double* g = new double[nmax];
	double* wa = new double[2 * mmax * nmax + 4 * nmax + 12 * mmax * mmax + 12 * mmax];


    //c     Declare a few additional variables for this sample problem.
    double t1, t2;
    int i;

    //c     We wish to have output at every iteration.
    iprint = -1;

    //c     We specify the tolerances in the stopping criteria.
    factr = 0;//1.e+1;
    pgtol = 0;//1.0e-9;

    //c     We specify the dimension n of the sample problem and the number
    //c        m of limited memory corrections stored.  (n and m should not
    //c        exceed the limits nmax and mmax respectively.)
    n = nmax;
    m = mmax;

    //c     We now provide nbd which defines the bounds on the variables:
    //c                    l   specifies the lower bounds,
    //c                    u   specifies the upper bounds. 

    //c     First set bounds on the odd-numbered variables. (*even-numbered in c*)
    for (i = 0; i < n; i += 2) {
        nbd[i] = 0;
        l[i] = 1.0;
        u[i] = 100.0;
    }

    //c     Next set bounds on the even-numbered variables. (*odd-numbered in c*)
    for (i = 1; i < n; i += 2) {
        nbd[i] = 0;
        l[i] = -100.0;
        u[i] = 100.0;
    }

    //c     We now define the starting point.
    for (i = 0; i < n; i++)
        x[i] = 3.0;

    printf("Solving sample problem.\n");
    printf(" (f=0.0 at the optimal solution.)\n");


    //c     We start the iteration by initializing task.
    // (**MUST clear remaining chars in task with spaces (else crash)!**)
    strcpy(task, "START");
    for (i = 5; i < SIXTY; i++)
        task[i] = ' ';

    //c     This is the call to the L-BFGS-B code.
    // (* call the L-BFGS-B routine with task='START' once before loop *)
    setulb_(&n, &m, x, l, u, nbd, &f, g, &factr, &pgtol, wa, iwa, task, &iprint,
            csave, lsave, isave, dsave, (long int) 60, (long int) 60);

    // (* while routine returns "FG" or "NEW_X" in task, keep calling it *)
    while (strncmp(task, "FG", 2) == 0 || strncmp(task, "NEW_X", 5) == 0) {

        if (strncmp(task, "FG", 2) == 0) {
            //c   the minimization routine has returned to request the
            //c   function f and gradient g values at the current x

            //c        Compute function value f for the sample problem.

            f = 0.25 * (x[0] - 1.0)*(x[0] - 1.0);
            for (i = 1; i < n; i++)
                f += pow(x[i] - x[i - 1] * x[i - 1], 2);
            f *= 4.0;

            //c        Compute gradient g for the sample problem.
            //t1 = (x[1]-x[0])*(x[1]-x[0]);
            t1 = x[1] - x[0] * x[0];
            g[0] = 2.0 * (x[0] - 1.0) - 16.0 * x[0] * t1;
            for (i = 1; i < n - 1; i++) {
                t2 = t1;
                //t1=(x[i+1]-x[i])*(x[i+1]-x[i]);
                t1 = x[i + 1] - x[i] * x[i];
                g[i] = 8.0 * t2 - 16.0 * x[i] * t1;
            }
            g[n - 1] = 8.0 * t1;

        } else {
            // the minimization routine has returned with a new iterate,
            // and we have opted to continue the iteration (do nothing here)
        }

        //c          go back to the minimization routine.
        setulb_(&n, &m, x, l, u, nbd, &f, g, &factr, &pgtol, wa, iwa, task, &iprint,
                csave, lsave, isave, dsave, (long int) 60, (long int) 60);
        cout<<f<<endl;
    }

    for (i = 0; i < 25; i++) {
        printf("i=%2d     x[i]=%f\n", i, x[i]);
        
    }
    
    //c     If task is neither FG nor NEW_X we terminate execution.
    // (* the minimization routine has returned with one of
    //  {CONV, ABNO, ERROR} *)

    return 0;
}